{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e87181a4-e19a-4b63-954a-64d4dadd9e22",
   "metadata": {},
   "source": [
    "# 载入code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1c415039-f322-47db-88d1-2aade69b5548",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Transcriptomic and Spatial Proteomic Profiling Reveals the Cellular Composition and Spatial Organization of the Human Bone Marrow Microenvironment\n",
    "# Author: Shovik Bandyopadhyay\n",
    "# This script documents the analysis of scRNA-Seq data in this paper for Figure 1,2 and 5F\n",
    "\n",
    "# Load necessary libraries ----\n",
    "# if (!require(\"BiocManager\", quietly = TRUE))\n",
    "#   install.packages(\"BiocManager\")\n",
    "# BiocManager::install(c(\"cowplot\"), force = TRUE)\n",
    "\n",
    "library(dplyr)\n",
    "library(Seurat)\n",
    "library(patchwork)  # 安装\n",
    "library(readr)\n",
    "library(RColorBrewer)\n",
    "library(ggplot2)\n",
    "library(VisCello)\n",
    "library(viridis)\n",
    "library(forcats)\n",
    "library(ggrastr)\n",
    "library(cowplot)\n",
    "library(GSEABase)\n",
    "library(SeuratDisk)\n",
    "library(ComplexHeatmap)\n",
    "library(DoubletFinder)\n",
    "library(irlba)\n",
    "library(AUCell)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59c45a77-cddf-4b85-9ccd-0fb0b1d96a74",
   "metadata": {},
   "source": [
    "# 双细胞处理函数\n",
    "\n",
    "* 注意v3 版本是旧版本，考虑参考文献更替成新版本"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60061df6-b589-418d-93ea-ed7efea246c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Doublet Finder ------\n",
    "FindDoublets <- function(seurat.rna, PCs = 1:50, exp_rate = 0.02, sct = FALSE){\n",
    "  # sct--do SCTransform or not\n",
    "  \n",
    "  ## pK identification\n",
    "  sweep.res.list <- paramSweep_v3(seurat.rna, PCs = PCs, sct = sct)\n",
    "  sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)\n",
    "  bcmvn <- find.pK(sweep.stats)\n",
    "  \n",
    "  ## Homotypic Doublet proportion Estimate\n",
    "  annotations <- seurat.rna@meta.data$seurat_clusters\n",
    "  homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- seu_kidney@meta.data$ClusteringResults\n",
    "  nExp_poi <- round(exp_rate * length(seurat.rna$seurat_clusters))  ## Assuming 7.5% doublet formation rate - tailor for your dataset\n",
    "  nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))\n",
    "  \n",
    "  ## Run DoubletFinder with varying classification stringencies ----------------------------------------------------------------\n",
    "  seurat.rna <- doubletFinder_v3(seurat.rna, PCs = PCs, pN = 0.25,\n",
    "                                 pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE, \n",
    "                                 sct = sct)\n",
    "  \n",
    "  seurat.rna <- doubletFinder_v3(seurat.rna, PCs = PCs, pN = 0.25, \n",
    "                                 pK = 0.09, nExp = nExp_poi.adj,\n",
    "                                 reuse.pANN = paste0(\"pANN_0.25_0.09_\", nExp_poi), \n",
    "                                 sct = sct)\n",
    "  doublet_var = paste0('DF.classifications_0.25_0.09_', nExp_poi.adj)\n",
    "  seurat.rna[['Doublet_Singlet']] = seurat.rna[[doublet_var]]\n",
    "  \n",
    "  mnames = names(seurat.rna@meta.data)\n",
    "  seurat.rna@meta.data[, grep(mnames, pattern = '0.25_0.09')] <- NULL\n",
    "  #seurat.rna = subset(seurat.rna, Doublet_Singlet == 'Singlet')\n",
    "  return(seurat.rna)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7db2c5c8-d53f-4b59-80fd-30a0bbefb44d",
   "metadata": {},
   "source": [
    "# Remove doublets, preprocess the data, compute dimension reductions and cluster ----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ebfdce-4aad-4ba5-8d9f-42170027baca",
   "metadata": {},
   "source": [
    "## Load the cellranger output files\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e652a741-63a6-4f84-a071-cf5cd1bdecb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "MACS <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB51_H21andH23_scRNA/H14/With_Introns/MACS/filtered_feature_bc_matrix/\")\n",
    "H21 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB51_H21andH23_scRNA/H21/filtered_feature_bc_matrix/\")\n",
    "H23 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB51_H21andH23_scRNA/H23/filtered_feature_bc_matrix/\")\n",
    "H24 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB62_scRNASeq_S2/H24/filtered_feature_bc_matrix/\")\n",
    "H32 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB62_scRNASeq_S2/H32/filtered_feature_bc_matrix/\")\n",
    "H33 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB62_scRNASeq_S2/H33/filtered_feature_bc_matrix/\")\n",
    "H34 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB62_scRNASeq_S2/H34/filtered_feature_bc_matrix/\")\n",
    "H35 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/SB62_scRNASeq_S2/H35/filtered_feature_bc_matrix/\")\n",
    "H36 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/NBM_Atlas_scRNA/H36/filtered_feature_bc_matrix/\")\n",
    "H38 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/NBM_Atlas_scRNA/H38/filtered_feature_bc_matrix/\")\n",
    "H39 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/NBM_Atlas_scRNA/H39/filtered_feature_bc_matrix/\")\n",
    "H41 <- Read10X(data.dir = \"~/Documents/NBM_Microenvironment/NBM_Atlas_scRNA/H41/filtered_feature_bc_matrix/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bac202f4-2ba2-4584-9ad4-8a1dc3632ca0",
   "metadata": {},
   "source": [
    "## Initialize the Seurat object with the raw (non-normalized data)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fe41ddd-f5da-4b88-b68d-7e66dfbd289a",
   "metadata": {},
   "outputs": [],
   "source": [
    "MACS <- CreateSeuratObject(counts = MACS, project = \"H14_MACS\", min.cells = 3, min.features = 100)\n",
    "MACS\n",
    "H21 <- CreateSeuratObject(counts = H21, project = \"H21\", min.cells = 3, min.features = 100)\n",
    "H21\n",
    "H23 <- CreateSeuratObject(counts = H23, project = \"H23\", min.cells = 3, min.features = 100)\n",
    "H23\n",
    "H24 <- CreateSeuratObject(counts = H24, project = \"H24\", min.cells = 3, min.features = 100)\n",
    "H24\n",
    "H32 <- CreateSeuratObject(counts = H32, project = \"H32\", min.cells = 3, min.features = 100)\n",
    "H32\n",
    "H33 <- CreateSeuratObject(counts = H33, project = \"H33\", min.cells = 3, min.features = 100)\n",
    "H33\n",
    "H34 <- CreateSeuratObject(counts = H34, project = \"H34\", min.cells = 3, min.features = 100)\n",
    "H34\n",
    "H35 <- CreateSeuratObject(counts = H35, project = \"H35\", min.cells = 3, min.features = 100)\n",
    "H35\n",
    "H36 <- CreateSeuratObject(counts = H36, project = \"H36\", min.cells = 3, min.features = 100)\n",
    "H36\n",
    "H38 <- CreateSeuratObject(counts = H38, project = \"H38\", min.cells = 3, min.features = 100)\n",
    "H38\n",
    "H39 <- CreateSeuratObject(counts = H39, project = \"H39\", min.cells = 3, min.features = 100)\n",
    "H39\n",
    "H41 <- CreateSeuratObject(counts = H41, project = \"H41\", min.cells = 3, min.features = 100)\n",
    "H41"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcee1838-f142-43fa-ab2e-bfd78e57d40a",
   "metadata": {},
   "source": [
    "## Preprocess samples individually to remove doublets -----   \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45014c4b-6f31-4ecf-9f44-2c5ab025213a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ob.list <- list(MACS, H21, H23,H24,H32,H33,H34,H35, H36, H38, H39, H41)\n",
    "for (i in 1:length(ob.list)) {\n",
    "  ob.list[[i]][[\"percent.mt\"]] <- PercentageFeatureSet(ob.list[[i]], pattern = \"^MT-\")\n",
    "  ob.list[[i]] <- subset(ob.list[[i]], subset = nFeature_RNA > 100 & nFeature_RNA < 10000 & percent.mt < 10)\n",
    "  # ob.list[[i]] <- subset(ob.list[[i]], cells = cells.use)\n",
    "  ob.list[[i]] <- NormalizeData(ob.list[[i]])\n",
    "  ob.list[[i]] <- FindVariableFeatures(ob.list[[i]])\n",
    "  ob.list[[i]] <- ScaleData(ob.list[[i]])\n",
    "  ob.list[[i]] <- RunPCA(ob.list[[i]], features = VariableFeatures(object = ob.list[[i]]))\n",
    "  ob.list[[i]] <- RunUMAP(ob.list[[i]], dims = 1:30, reduction.name = \"UMAP_dim30\", reduction.key = \"UMAP_dim30_\")\n",
    "  ob.list[[i]] <- RunUMAP(ob.list[[i]], dims = 1:50, reduction.name = \"UMAP_dim50\", reduction.key = \"UMAP_dim50_\")\n",
    "  ob.list[[i]] <- FindNeighbors(ob.list[[i]], dims = 1:30)\n",
    "  ob.list[[i]] <- FindClusters(ob.list[[i]], resolution = 1)\n",
    "  ob.list[[i]] <- FindDoublets(ob.list[[i]], PCs = 1:30, sct = FALSE, exp_rate = (length(colnames(ob.list[[i]]))/125000))\n",
    "}\n",
    "\n",
    "MACS <- ob.list[[1]]\n",
    "H21 <- ob.list[[2]]\n",
    "H23 <- ob.list[[3]]\n",
    "H24 <- ob.list[[4]]\n",
    "H32 <- ob.list[[5]]\n",
    "H33 <- ob.list[[6]]\n",
    "H34 <- ob.list[[7]]\n",
    "H35 <- ob.list[[8]]\n",
    "H36 <- ob.list[[9]]\n",
    "H38 <- ob.list[[10]]\n",
    "H39 <- ob.list[[11]]\n",
    "H41 <- ob.list[[12]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b78a8507-395a-47b2-a7d4-e8babb022eef",
   "metadata": {},
   "source": [
    "## Merge all seurat objects together "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "587bb617-98d6-4f93-92ed-58db49c12615",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined <- merge(x = MACS, y = c(H21, H23,H24, H32, H33,H34,H35,H36, H38, H39, H41), add.cell.ids = c(\"H14_MACS\", \"H21\", \"H23\",\"H24\", \"H32\", \"H33\", \"H34\", \"H35\", \"H36\", \"H38\", \"H39\",\"H41\"), project = \"SB62_NBM\")\n",
    "\n",
    "ob.list <- list(combined)\n",
    "\n",
    "for (i in 1:length(ob.list)) {\n",
    "  ob.list[[i]][[\"percent.mt\"]] <- PercentageFeatureSet(ob.list[[i]], pattern = \"^MT-\")\n",
    "  ob.list[[i]] <- subset(ob.list[[i]], subset = nFeature_RNA > 100 & nFeature_RNA < 10000 & percent.mt < 10)\n",
    "  ob.list[[i]] <- subset(ob.list[[i]], subset = Doublet_Singlet == \"Singlet\")\n",
    "  ob.list[[i]] <- NormalizeData(ob.list[[i]])\n",
    "  ob.list[[i]] <- FindVariableFeatures(ob.list[[i]])\n",
    "  ob.list[[i]] <- ScaleData(ob.list[[i]])\n",
    "  ob.list[[i]] <- RunPCA(ob.list[[i]], features = VariableFeatures(object = ob.list[[i]]))\n",
    "  ob.list[[i]] <- RunUMAP(ob.list[[i]], dims = 1:30, reduction.name = \"UMAP_dim30\", reduction.key = \"UMAP_dim30_\")\n",
    "  ob.list[[i]] <- RunUMAP(ob.list[[i]], dims = 1:50, reduction.name = \"UMAP_dim50\", reduction.key = \"UMAP_dim50_\")\n",
    "  ob.list[[i]] <- FindNeighbors(ob.list[[i]], dims = 1:30)\n",
    "  ob.list[[i]] <- FindClusters(ob.list[[i]], algorithm = 2, resolution = 1)\n",
    "}\n",
    "\n",
    "combined <- FindClusters(combined, algorithm = 2, resolution = 1.5)\n",
    "\n",
    "combined <- ob.list[[1]] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c57d83c-09ca-4990-bf5e-15d4617042a1",
   "metadata": {},
   "source": [
    "## remove contaminating cells and doublets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a78caff2-aac2-402a-b8f9-297c264048db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find and save all markers for combined UMAP\n",
    "combined_markers <- FindAllMarkers(combined, max.cells.per.ident = 1000) # Downsample 1000 cells per cluster\n",
    "write.csv(combined_markers, \"SB66_combined_markers.csv\")\n",
    "\n",
    "# Cluster 46 - muscle\n",
    "# Cluster 41 - CXCL14+ Fibroblasts\n",
    "# Cluster 45 - Doublets Plasma + T\n",
    "combined <- readRDS(\"SB66_combined.RDS\")\n",
    "combined <- subset(combined, subset = seurat_clusters != 46 & seurat_clusters != 41 & seurat_clusters != 45)\n",
    "\n",
    "ob.list <- list(combined)\n",
    "\n",
    "for (i in 1:length(ob.list)) {\n",
    "  ob.list[[i]][[\"percent.mt\"]] <- PercentageFeatureSet(ob.list[[i]], pattern = \"^MT-\")\n",
    "  ob.list[[i]] <- subset(ob.list[[i]], nFeature_RNA > 300 & nCount_RNA > 1000 & nFeature_RNA < 10000 & percent.mt < 10)\n",
    "  ob.list[[i]] <- subset(ob.list[[i]], subset = Doublet_Singlet == \"Singlet\")\n",
    "  ob.list[[i]] <- NormalizeData(ob.list[[i]])\n",
    "  ob.list[[i]] <- FindVariableFeatures(ob.list[[i]])\n",
    "  ob.list[[i]] <- ScaleData(ob.list[[i]])\n",
    "  ob.list[[i]] <- RunPCA(ob.list[[i]], features = VariableFeatures(object = ob.list[[i]]))\n",
    "  ob.list[[i]] <- RunUMAP(ob.list[[i]], dims = 1:30, reduction.name = \"UMAP_dim30\", reduction.key = \"UMAP_dim30_\")\n",
    "  ob.list[[i]] <- RunUMAP(ob.list[[i]], dims = 1:50, reduction.name = \"UMAP_dim50\", reduction.key = \"UMAP_dim50_\")\n",
    "  ob.list[[i]] <- FindNeighbors(ob.list[[i]], dims = 1:30)\n",
    "  ob.list[[i]] <- FindClusters(ob.list[[i]], algorithm = 2, resolution = 1.5)\n",
    "}\n",
    "\n",
    "\n",
    "combined <- ob.list[[1]] \n",
    "\n",
    "# Find and save all markers for combined UMAP\n",
    "combined_markers <- FindAllMarkers(combined, max.cells.per.ident = 1000) # Downsample 1000 cells per cluster\n",
    "write.csv(combined_markers, \"SB66_combined_markers.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cebb8653-f931-47c0-82c9-a74d89c5b2e7",
   "metadata": {},
   "source": [
    "# Annotate the cells in the final dataset ----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fa29701-1af5-448a-9b60-f584ceb1264a",
   "metadata": {},
   "source": [
    "## Prepare seurat object for viscello"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "740f0aaf-bba5-463a-b2e0-fa505f8fb622",
   "metadata": {},
   "outputs": [],
   "source": [
    "library(VisCello)\n",
    "DefaultAssay(combined) <- \"RNA\"\n",
    "fmeta <- data.frame(symbol = rownames(combined)) \n",
    "rownames(fmeta) <- fmeta$symbol\n",
    "eset <- new(\"ExpressionSet\",\n",
    "            assayData = assayDataNew(\"environment\", exprs=combined@assays$RNA@counts, \n",
    "                                     norm_exprs = combined@assays$RNA@data),\n",
    "            phenoData =  new(\"AnnotatedDataFrame\", data = combined@meta.data),\n",
    "            featureData = new(\"AnnotatedDataFrame\", data = fmeta))\n",
    "saveRDS(eset, \"VisCello/eset.rds\") \n",
    "# Creating a cello for all the cells\n",
    "cello <- new(\"Cello\", name = \"Combined All Cells SB66\", idx = 1:ncol(eset)) # Index is basically the column index, here all cells are included \n",
    "# Code for computing dimension reduction, not all of them is necessary, and you can input your own dimension reduction result into the cello@proj list.\n",
    "# It is also recommended that you first filter your matrix to remove low expression genes and cells, and input a matrix with variably expressed genes\n",
    "cello@proj <- list('PCA' = combined@reductions$pca@cell.embeddings,\n",
    "                   'UMAP_dim30' =combined@reductions$UMAP_dim30@cell.embeddings,\n",
    "                   'UMAP_dim50' =combined@reductions$UMAP_dim50@cell.embeddings)\n",
    "clist <- list()\n",
    "clist[[\"SB66 Global dataset\"]] <- cello\n",
    "saveRDS(clist, \"VisCello/clist.rds\") \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21a99bea-b785-45eb-9765-30af9d7e4481",
   "metadata": {},
   "source": [
    "## Rename clusters after manual inspection of viscello obj and markers csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45cb28de-dd7d-47d0-89a5-023335830783",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined <- SetIdent(combined, value = \"seurat_clusters\")\n",
    "new.cluster.ids.combined <- c(\"Adipo-MSC\", \"Plasma Cell\", \"Plasma Cell\", \"Plasma Cell\", \"THY1+ MSC\", \"Osteo-MSC\",\n",
    "                              \"Neutrophil\", \"SEC\", \"Pro-B\", \"CD4+ T-Cell\", \"Osteoblast\", \"VSMC\", \"HSC\",\n",
    "                              \"Late Myeloid\", \"Late Erythroid\", \"GMP\", \"MEP\", \"OsteoFibro-MSC\",\n",
    "                              \"CD8+ T-Cell\", \"MPP\", \"Early Myeloid Progenitor\", \"GMP\", \"Fibro-MSC\", \"RNAlo MSC\", \"Monocyte\", \"Mature B\",\n",
    "                              \"AEC\", \"Cycling HSPC\", \"Pre-Pro B\", \"Late Myeloid\", \"Erythroblast\", \"Plasma Cell\", \"VSMC\",\n",
    "                              \"Pre-B\", \"pDC\", \"Ba/Eo/Ma\", \"AEC\", \"CLP\", \"Plasma Cell\",\n",
    "                              \"Cycling DCs\", \"RBC\", \"Plasma Cell\", \"MPP\", \"Macrophages\", \"RBC\", \"NKT Cell\")\n",
    "\n",
    "names(new.cluster.ids.combined) <- levels(combined)\n",
    "combined <- RenameIdents(combined, new.cluster.ids.combined)\n",
    "DimPlot(combined, label = TRUE, reduction = \"UMAP_dim30\") + coord_fixed() + NoAxes() + NoLegend()\n",
    "combined$cluster_anno_l2 <- combined@active.ident"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1ca6951-a692-4a11-8a1b-abe41407fa76",
   "metadata": {},
   "source": [
    "## Mks do not cluster separately, need to subcluster erythroblasts to properly annotate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "433df73b-51b7-4b74-a4f0-fcb6a10ddb27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset erythroblast cluster and rename Mks \n",
    "EBs <- subset(combined, cluster_anno_l2 == \"Erythroblast\")\n",
    "EBs <- NormalizeData(EBs)\n",
    "EBs <- FindVariableFeatures(EBs)\n",
    "EBs <- ScaleData(EBs)\n",
    "EBs <- RunPCA(EBs, reduction.name = \"EB_pca\", features = VariableFeatures(EBs))\n",
    "EBs <- RunUMAP(EBs, dims = 1:10, reduction = \"EB_pca\", reduction.name = \"EB_UMAP_dim10\")\n",
    "EBs <- RunUMAP(EBs, dims = 1:20, reduction = \"EB_pca\", reduction.name = \"EB_UMAP_dim20\")\n",
    "EBs <- RunUMAP(EBs, dims = 1:30, reduction = \"EB_pca\", reduction.name = \"EB_UMAP_dim30\")\n",
    "EBs <- RunUMAP(EBs, dims = 1:40, reduction = \"EB_pca\", reduction.name = \"EB_UMAP_dim40\")\n",
    "EBs <- RunUMAP(EBs, dims = 1:50, reduction = \"EB_pca\", reduction.name = \"EB_UMAP_dim50\")\n",
    "EBs <- FindNeighbors(EBs, dims = 1:20)\n",
    "EBs <- FindClusters(EBs, algorithm = 2, resolution = 0.5, group.singletons = FALSE)\n",
    "DimPlot(EBs, label = TRUE, reduction = \"EB_UMAP_dim50\") + coord_fixed() + NoAxes() + NoLegend()\n",
    "FeaturePlot(EBs, features = c(\"ITGA2B\", \"GATA2\"), reduction = \"EB_UMAP_dim50\", coord.fixed = TRUE) & NoAxes() \n",
    "VlnPlot(EBs, features = c(\"ITGA2B\",\"ITGB3\", \"GATA2\" ,\"GYPC\")) # can see that cluster 4 is most megakaryocytic\n",
    "# Rename Mk cluster based on CD41 expression\n",
    "Mk_ids <- colnames(subset(EBs, subset = seurat_clusters == 4)) # GATA2/CD41/CD61 expression\n",
    "combined$cluster_anno_l2 <- as.character(combined$cluster_anno_l2)\n",
    "combined$cluster_anno_l2[Mk_ids] <- \"Megakaryocyte\" \n",
    "combined$cluster_anno_l2 <- as.factor(combined$cluster_anno_l2)\n",
    "levels(combined$cluster_anno_l2)\n",
    "DimPlot(combined, label= TRUE, repel = TRUE, group.by = 'cluster_anno_l2') + NoAxes() + coord_fixed()\n",
    "\n",
    "# remove NKT cluster because they are likely doublets (T-cell markers, Osteolineage markers co-expressed)\n",
    "combined <- subset(combined, cluster_anno_l2 != \"NKT Cell\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35cab787-52a3-4db8-87e1-c3c35da17be2",
   "metadata": {},
   "source": [
    "## Make coarse annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9206c57c-2bd4-4535-9a26-719210f25db5",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cluster_anno_coarse <- c(\"Mesenchymal\", \"Endothelial\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Mesenchymal\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Hematopoietic\", \"Mesenchymal\", \"Endothelial\", \"Mesenchymal\", \"Muscle\")\n",
    "combined <- SetIdent(combined, value = \"cluster_anno_l2\")\n",
    "names(cluster_anno_coarse) <- levels(combined)\n",
    "combined <- RenameIdents(combined, cluster_anno_coarse)\n",
    "DimPlot(combined, label = TRUE, reduction = \"UMAP_dim30\") + coord_fixed() + NoAxes() + NoLegend()\n",
    "combined$cluster_anno_coarse <- combined@active.ident\n",
    "coarse_cols <- c('#f88379','#FFA750','#6495ED','#B8242D')\n",
    "DimPlot(combined, label = TRUE, reduction = \"UMAP_dim30\", cols=coarse_cols) + coord_fixed() + NoAxes() + NoLegend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4b4a473-d5ac-402b-8d9b-1da64caf5fde",
   "metadata": {},
   "source": [
    "## save cluster_anno_l2 colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b2e57fc-1dd1-4d73-ac91-a2760d3e9628",
   "metadata": {},
   "outputs": [],
   "source": [
    "cal2_cols <- c(\"#FFD580\", \"#FFBF00\", \"#B0DFE6\", \"#7DB954\", \"#64864A\", \"#8FCACA\", \"#4682B4\", \"#CAA7DD\", \"#B6D0E2\", \"#a15891\", \"#FED7C3\", \"#A8A2D2\", \"#CF9FFF\", \"#9C58A1\", \"#2874A6\", \"#96C5D7\", \"#63BA97\", \"#BF40BF\", \"#953553\", \"#6495ED\", \"#E7C7DC\", \"#5599C8\", \"#FA8072\", \"#F3B0C3\", \"#F89880\", \"#40B5AD\", \"#019477\", \"#97C1A9\", \"#C6DBDA\", \"#CCE2CB\", \"#79127F\", \"#FFC5BF\", \"#ff9d5c\", \"#FFC8A2\", \"#DD3F4E\")\n",
    "cal2_col_names <- c(\"Adipo-MSC\", \"AEC\", \"Ba/Eo/Ma\", \"CD4+ T-Cell\", \"CD8+ T-Cell\", \"CLP\", \"Cycling DCs\", \"Cycling HSPC\", \"Early Myeloid Progenitor\", \"Erythroblast\", \"Fibro-MSC\", \"GMP\", \"HSC\", \"Late Erythroid\", \"Late Myeloid\", \"Macrophages\", \"Mature B\", \"Megakaryocyte\", \"MEP\", \"Monocyte\", \"MPP\", \"Neutrophil\", \"Osteo-MSC\", \"Osteoblast\", \"OsteoFibro-MSC\", \"pDC\", \"Plasma Cell\", \"Pre-B\", \"Pre-Pro B\", \"Pro-B\", \"RBC\", \"RNAlo MSC\", \"SEC\", \"THY1+ MSC\", \"VSMC\")\n",
    "names(cal2_cols) <- cal2_col_names\n",
    "\n",
    "combined <- SetIdent(combined, value = \"cluster_anno_l2\")\n",
    "combined_markers_anno <- FindAllMarkers(combined, max.cells.per.ident = 500) # Downsample 1000 cells per cluster\n",
    "write.csv(combined_markers_anno, \"SB66_combined_markers_annotated.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e0ae34c-1b89-4918-8102-210b30cda5a9",
   "metadata": {},
   "source": [
    "# Figure 1 Generation ----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba8a7b12-2a20-4885-9ddb-bf0470edd49a",
   "metadata": {},
   "source": [
    "## in text QC metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "014c5f51-ea53-46c2-962f-d063afe8ebc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined@meta.data %>% dplyr::summarise(c(mean(nCount_RNA), mean(nFeature_RNA), mean(percent.mt)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b12b0ec-5cd6-4aee-8e4f-4786dde80ee7",
   "metadata": {},
   "source": [
    "## Supplemental Figure S1C QC Metrics -----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd60bf26-1333-4a84-a25b-9625fba37c88",
   "metadata": {},
   "outputs": [],
   "source": [
    "VlnPlot(combined, group.by =\"cluster_anno_coarse\", features = c(\"nFeature_RNA\", \"nCount_RNA\", \"percent.mt\"), pt.size=0) & scale_fill_manual(values =c(\"#AEC7E8\", \"#FFBB78\", \"#98DF8A\", \"#FF9896\"))\n",
    "VlnPlot(combined, group.by =\"cluster_anno_coarse\", features = c(\"nCount_RNA\"), pt.size=0, y.max = quantile(combined$nCount_RNA,0.99)) & scale_fill_manual(values =c(\"#AEC7E8\", \"#FFBB78\", \"#98DF8A\", \"#FF9896\"))\n",
    "\n",
    "combined@meta.data %>% dplyr::group_by(cluster_anno_coarse)  %>% summarise(median(nFeature_RNA), median(nCount_RNA), median(percent.mt)) %>% gt() -> supp_fig_1D\n",
    "gtsave(supp_fig_1D, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure1/Supp_Fig_1D.pdf\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "963bf648-9424-4dd7-b787-490bf742efed",
   "metadata": {},
   "source": [
    "## Supplemental Figure S1D UMAP Feature Plots -----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0debb782-d6cf-4210-a745-f0201a4c9d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 <- FeaturePlot(object = combined, raster = TRUE, raster.dpi = c(1028,1028), features = c(\"CXCL12\",\"NCAM1\",\"CDH5\", \"PTPRC\", \"MZB1\", \"CSF3R\"), cols = brewer.pal(n = 100, name = \"Reds\"),ncol = 6, max.cutoff = 'q99', coord.fixed = TRUE) & NoAxes() \n",
    "p1_raster <- rasterize(p1)\n",
    "ggsave(p1_raster,file = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure1/PanelE_scRNA_MarkerFeaturePlots.pdf\", device = \"pdf\", width = 15, height = 2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efad5c6f-e0c4-4307-ad4b-dfeba50c0033",
   "metadata": {},
   "source": [
    "## Figure S1E CODEX Cell Type Frequencies Per Sample By Age ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "137c1303-2bc1-49fa-8e75-bd1d1fada6a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(data=percentage_data, aes(x=Age,y=percentage, fill=cluster_anno_l1)) +\n",
    "  geom_bar(position=\"fill\",stat=\"identity\") + \n",
    "  theme_bw() + \n",
    "  labs(y='Cell Type Frequency') + \n",
    "  labs(x='') + scale_fill_manual(values=cal1_cols) + theme(axis.text = element_text(size=16)) + ggtitle(\"CODEX Cell Type Frequencies Per Sample\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f87b3ac5-f9b9-451f-b72d-cb20aee533fa",
   "metadata": {},
   "source": [
    "## perform linear regression testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16dbba9d-302f-46a5-8a65-cef5f1f2d781",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "linear_regression_results <- percentage_data %>%\n",
    "  group_by(cluster_anno_l1) %>%\n",
    "  do(model = lm(percentage ~ Age, data = .))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69b5258d-d95d-4580-ac68-bb6b272105a3",
   "metadata": {},
   "source": [
    "## Extract p-values and R-squared values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ef2c8f-39e4-4198-ae4e-ba680d6e446b",
   "metadata": {},
   "outputs": [],
   "source": [
    "result_df <- linear_regression_results %>%\n",
    "  summarise(\n",
    "    Cell_Type = cluster_anno_l1,\n",
    "    P_Value = coef(summary(model))[, \"Pr(>|t|)\"][2],\n",
    "    R_Value = summary(model)$r.squared\n",
    "  )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44afdb98-dee6-43f7-86ba-a3368884eb87",
   "metadata": {},
   "source": [
    "## confirm everything adds up to 100 per patient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d06f0b13-df54-43be-b52b-b8ad5d106519",
   "metadata": {},
   "outputs": [],
   "source": [
    "percentage_data %>% group_by(orig.ident) %>% summarise(sum(percentage))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "029d8ba3-5e61-4424-95d1-ca76c9eeee09",
   "metadata": {},
   "source": [
    "# repeat for MSCs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84cd1c63-1c46-4703-b852-73710f465321",
   "metadata": {},
   "source": [
    "## plot cell types as function of age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d98b4e2f-3eb7-4568-9b2d-71b901b4e6fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# \n",
    "percentage_data <- MSCs@meta.data %>%\n",
    "  group_by(orig.ident, cluster_anno_l2) %>%\n",
    "  summarise(count = n()) %>% \n",
    "  mutate(percentage = count / sum(count) * 100)\n",
    "\n",
    "percentage_data <- left_join(percentage_data, additional_md, by = \"orig.ident\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4fd2828-1f62-41c0-8be5-0364599ada4b",
   "metadata": {},
   "source": [
    "## Create the ggplot2 plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f69c6f04-bb62-4ae0-94ff-2b0c75ebcc2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(percentage_data, aes(x = Age, y = percentage, color = cluster_anno_l2)) +\n",
    "  geom_line() + geom_point() +\n",
    "  labs(title = \"Cell Type Percentages by Age\",\n",
    "       x = \"Age\",\n",
    "       y = \"Percentage\") +\n",
    "  scale_fill_brewer(palette = \"Set1\") +  # Choose a color palette\n",
    "  theme_minimal()\n",
    "\n",
    "linear_regression_results <- percentage_data %>%\n",
    "  group_by(cluster_anno_l2) %>%\n",
    "  do(model = lm(percentage ~ Age, data = .))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce198fe-6edb-4621-8c5b-9055c1a0d0d2",
   "metadata": {},
   "source": [
    "## Extract p-values and R-squared values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8360f9d6-2170-48fd-bc7e-bac2ffae66af",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "result_df <- linear_regression_results %>%\n",
    "  summarise(\n",
    "    Cell_Type = cluster_anno_l2,\n",
    "    P_Value = coef(summary(model))[, \"Pr(>|t|)\"][2],\n",
    "    R_Value = summary(model)$r.squared\n",
    "  )\n",
    "print(result_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a63305e-f562-4130-ac42-4ccfc31d7b82",
   "metadata": {},
   "source": [
    "## Supplemental Figure S1F Azimuth Comparison ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad4b212a-bd18-42c5-b452-1102cc686ada",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "azimuth <- readRDS(\"~/Documents/AML_Microenvironment/SB36_CibersortX/Azimuth_BM_Reference/ref.Rds\")\n",
    "DimPlot(azimuth, reduction = \"refUMAP\", group.by = \"celltype.l2\") + coord_fixed() + NoAxes()\n",
    "table(azimuth$celltype.l2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63a34902-dabd-480c-bef5-aac2d4dd7f07",
   "metadata": {},
   "source": [
    "# Figure 1B - UMAP Showing all of the Cell Types Captured in Our Analysis ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "203f3ba0-fe7f-48f6-8193-518dcc7298ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 <- DimPlot(combined,group.by = 'cluster_anno_l2',raster = TRUE, raster.dpi = c(1028,1028),label = TRUE, repel = TRUE, reduction = \"UMAP_dim30\", cols=cal2_cols) + coord_fixed() + NoAxes() + NoLegend()\n",
    "p1_raster <- rasterize(p1)\n",
    "ggsave(p1_raster,file = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Figure1/PanelB_scRNA_UMAP.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61f68489-417d-401e-8f44-dadbed484987",
   "metadata": {},
   "source": [
    "# Change levels "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6025d106-3449-4d67-b43c-394f1602a220",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cell_order <- c(\"HSC\", \"MPP\",\"Cycling HSPC\", \"MEP\", \"Erythroblast\", \"Late Erythroid\", \"RBC\",\"Megakaryocyte\",\"GMP\", \"Early Myeloid Progenitor\", \"Late Myeloid\", \"Neutrophil\", \"Monocyte\",\"Macrophages\",\"Ba/Eo/Ma\", \"Cycling DCs\", \"pDC\", \"CLP\", \"Pre-Pro B\", \"Pro-B\", \"Pre-B\", \"Mature B\", \"Plasma Cell\",\"CD4+ T-Cell\", \"CD8+ T-Cell\", \"NKT Cell\", \"AEC\", \"SEC\", \"VSMC\", \"THY1+ MSC\", \"Adipo-MSC\", \"OsteoFibro-MSC\", \"Osteo-MSC\", \"Osteoblast\", \"Fibro-MSC\", \"RNAlo MSC\")\n",
    "combined$cluster_anno_l2 <- factor(combined$cluster_anno_l2, levels = cell_order)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bed6b057-edba-427b-8bfc-e38dafdfdeba",
   "metadata": {},
   "source": [
    "# Figure 1C Create barplot showing lineage frequencies ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a3e4cbd-bf6d-4e92-a567-41c1e3cb20cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "## Create the level 1 (coarse) annotation based on broad lineages\n",
    "# cluster_anno_l1 <- c(\"Mesenchymal\", \"Endothelial\", \"Myeloid\", \"Lymphoid\", \"Lymphoid\", \"HSPC\",\"Myeloid\",\"HSPC\", \"Myeloid\", \"Meg/E\", \"Mesenchymal\",\"HSPC\",\"HSPC\",\"Meg/E\",\"Myeloid\",\"Myeloid\",\"Lymphoid\",\"Meg/E\",\"HSPC\",\"Myeloid\",\"HSPC\",\"Myeloid\",\"Mesenchymal\",\"Mesenchymal\",\"Mesenchymal\",\"Myeloid\",\"Lymphoid\",\"Lymphoid\",\"Lymphoid\",\"Lymphoid\",\"Meg/E\", \"Mesenchymal\", \"Endothelial\",  \"Mesenchymal\", \"Muscle\") # First make level 1 annotation resolved for the major hematopoietic lineages\n",
    "cluster_anno_l1 <- c(\"HSPC\", \"HSPC\", \"HSPC\", \"Meg/E\", \"Meg/E\", \"Meg/E\", \"Meg/E\", \"Meg/E\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Endothelial\", \"Endothelial\", \"Muscle\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\")\n",
    "combined <- SetIdent(combined, value = \"cluster_anno_l2\")\n",
    "names(cluster_anno_l1) <- levels(combined)\n",
    "combined <- RenameIdents(combined, cluster_anno_l1)\n",
    "combined$cluster_anno_l1 <- combined@active.ident\n",
    "combined$cluster_anno_l1 <- factor(combined$cluster_anno_l1, levels = c(\"HSPC\", \"Myeloid\", \"Lymphoid\", \"Meg/E\", \"Mesenchymal\", \"Endothelial\", \"Muscle\"))\n",
    "cal1_cols <- c(\"#E0B0FF\", \"#A7C7E7\", \"#AFE1AF\", \"#BDB5D5\", \"#FFB6C1\", \"#F28C28\", \"#DD3F4E\")\n",
    "DimPlot(combined, group.by = 'cluster_anno_l1',cols = cal1_cols, label = TRUE, reduction = \"UMAP_dim30\") + coord_fixed() + NoAxes() + NoLegend() # Just make sure reannotation has been done properly\n",
    "\n",
    "combined$cluster_anno_l2 <- droplevels(combined$cluster_anno_l2)\n",
    "col_coarse <- c(\"#FFB6C1\", \"#FFBF00\", \"#A7C7E7\", \"#AFE1AF\", \"#AFE1AF\", \"#AFE1AF\", \"#40B5AD\", \"#E0B0FF\", \"#A7C7E7\", \"#BDB5D5\", \"#FFB6C1\", \"#A7C7E7\", \"#E0B0FF\", \"#BDB5D5\", \"#A7C7E7\", \"#96C5D7\", \"#AFE1AF\", \"#CAA7DD\", \"#BDB5D5\", \"#6495ED\", \"#E0B0FF\", \"#A7C7E7\", \"#c3cede\", \"#c3cede\", \"#FFB6C1\", \"#40B5AD\", \"#AFE1AF\", \"#AFE1AF\", \"#AFE1AF\", \"#AFE1AF\", \"#BDB5D5\", \"#FFB6C1\", \"#ff9d5c\", \"#FFB6C1\", \"#DD3F4E\")\n",
    "names_coarse <- c(\"Mesenchymal\", \"AEC\", \"Myeloid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"pDC\", \"HSPC\", \"Myeloid\", \"Erythroid\", \"Mesenchymal\", \"Myeloid\", \"HSPC\", \"Erythroid\", \"Myeloid\", \"Macrophage\", \"Lymphoid\", \"Megakaryocyte\", \"Erythroid\", \"Monocyte\", \"HSPC\", \"Myeloid\", \"Osteolineage\", \"Osteolineage\", \"Mesenchymal\", \"pDC\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Erythroid\", \"Mesenchymal\", \"SEC\", \"Mesenchymal\", \"VSMC\")\n",
    "names(col_coarse) <- names_coarse\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19d3fa30-97fc-4435-bd76-dcf09dbbfba2",
   "metadata": {},
   "source": [
    "## Figure 1C - Bar chart showing the counts of each lineage assayed -----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "026708eb-2eba-4119-965f-e5b53be6cd9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Figure 1C (left)\n",
    "p2 <- combined@meta.data %>% ggplot(aes(y=forcats::fct_rev(forcats::fct_infreq(cluster_anno_l1)), fill = cluster_anno_l1)) + geom_bar(stat = 'count') +\n",
    "  theme_bw() + \n",
    "  scale_fill_manual(values = cal1_cols, labels = c(\"HSPC\", \"Myeloid\", \"Lymphoid\", \"Meg/E\", \"Mesenchymal\", \"Endothelial\", \"Muscle\")) + \n",
    "  theme(axis.text = element_text(size=12)) \n",
    "\n",
    "cell_counts <- as.data.frame(table(combined$cluster_anno_l1,combined$orig.ident))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4163bf26-ade1-4123-b526-4b3a015fa390",
   "metadata": {},
   "source": [
    "# Figure 1D Heatmap ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6ec0830-b71f-46ad-b982-f6006e88517b",
   "metadata": {},
   "outputs": [],
   "source": [
    "p3 <- ggplot(data=cell_counts, aes(x=forcats::fct_rev(Var2),y=Freq, fill=Var1)) +\n",
    "  geom_bar(position=\"fill\",stat=\"identity\") + \n",
    "  coord_flip() +\n",
    "  theme_bw() + \n",
    "  labs(y='Cell Type Frequency') + \n",
    "  labs(x='') + scale_fill_manual(values=cal1_cols) + theme(axis.text = element_text(size=16)) + ggtitle(\"CODEX Cell Type Frequencies Per Sample\") + scale_x_discrete(labels = c(\"H41\", \"H39\", \"H38\", \"H36\", \"H35\", \"H34\", \"H33\", \"H32\", \"H24\", \"H23\", \"H21\", \"H14\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1213fcdc-29e7-4472-8c9d-714992247813",
   "metadata": {},
   "source": [
    "## Figure 1D Heatmap ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04974236-0420-45e1-8511-5046f1c78fb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "library(ComplexHeatmap)\n",
    "library(circlize)\n",
    "combined_markers_anno %>%\n",
    "  group_by(cluster) %>%\n",
    "  top_n(n = 10, wt = 1/p_val_adj) -> top10 # Get top 10 most significant genes for each cell type\n",
    "combined <- SetIdent(combined, value = \"cluster_anno_l2\")\n",
    "avg <- AverageExpression(object = combined, group.by = 'cluster_anno_l2', slot = 'data',features = top10$gene) # Return average expression values across cells in each cluster for the selected genes from top3 object\n",
    "\n",
    "col_fun = colorRamp2(c(-4, 0, 4), c(\"blue\", \"white\", \"red\")) \n",
    "cell_order <- c(\"HSC\", \"MPP\",\"Cycling HSPC\", \"MEP\", \"Erythroblast\", \"Late Erythroid\", \"RBC\",\"Megakaryocyte\",\"GMP\", \"Early Myeloid Progenitor\", \"Late Myeloid\", \"Neutrophil\", \"Monocyte\",\"Macrophages\",\"Ba/Eo/Ma\", \"Cycling DCs\", \"pDC\", \"CLP\", \"Pre-Pro B\", \"Pro-B\", \"Pre-B\", \"Mature B\", \"Plasma Cell\",\"CD4+ T-Cell\", \"CD8+ T-Cell\", \"AEC\", \"SEC\", \"VSMC\", \"THY1+ MSC\", \"Adipo-MSC\", \"OsteoFibro-MSC\", \"Osteo-MSC\", \"Osteoblast\", \"Fibro-MSC\", \"RNAlo MSC\")\n",
    "cell_lineages <- c(\"HSPC\", \"HSPC\", \"HSPC\", \"Meg/E\", \"Meg/E\", \"Meg/E\", \"Meg/E\",\"Meg/E\",\"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Myeloid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Lymphoid\", \"Endothelial\", \"Endothelial\", \"Muscle\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\", \"Mesenchymal\",\"Mesenchymal\")\n",
    "cell_lineage_colors <- c(\"HSPC\" = \"steelblue\", \"HSPC\"= \"steelblue\", \"HSPC\"= \"steelblue\", \"Meg/E\" = 'darkorchid1', \"Meg/E\" = 'darkorchid1', \"Meg/E\" = 'darkorchid1', \"Meg/E\" = 'darkorchid1', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Myeloid\" = 'salmon', \"Lymphoid\" = 'darkgoldenrod1', \"Lymphoid\"= 'darkgoldenrod1', \"Lymphoid\"= 'darkgoldenrod1', \"Lymphoid\" = 'darkgoldenrod1', \"Lymphoid\" = 'darkgoldenrod1', \"Lymphoid\" = 'darkgoldenrod1', \"Lymphoid\"= 'darkgoldenrod1', \"Lymphoid\"= 'darkgoldenrod1', \"Endothelial\" = 'lavender', \"Endothelial\" = 'lavender', \"Muscle\" = 'cornflowerblue', \"Mesenchymal\" = 'darkolivegreen2', \"Mesenchymal\" = 'darkolivegreen2', \"Mesenchymal\" = 'darkolivegreen2', \"Mesenchymal\" = 'darkolivegreen2', \"Mesenchymal\" = 'darkolivegreen2', \"Mesenchymal\" = 'darkolivegreen2',\"Mesenchymal\" = 'darkolivegreen2')\n",
    "cell_lineage_colors <- c(\"HSPC\" = \"#E0B0FF\", \"HSPC\"= \"#E0B0FF\", \"HSPC\"= \"#E0B0FF\", \"Meg/E\" = \"#BDB5D5\", \"Meg/E\" = \"#BDB5D5\", \"Meg/E\" = \"#BDB5D5\", \"Meg/E\" = \"#BDB5D5\", \"Meg/E\" = \"#BDB5D5\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Myeloid\" = \"#A7C7E7\", \"Lymphoid\" = \"#AFE1AF\", \"Lymphoid\"= \"#AFE1AF\", \"Lymphoid\"= \"#AFE1AF\", \"Lymphoid\" = \"#AFE1AF\", \"Lymphoid\" = \"#AFE1AF\", \"Lymphoid\" = \"#AFE1AF\", \"Lymphoid\"= \"#AFE1AF\", \"Lymphoid\"= \"#AFE1AF\", \"Endothelial\" = \"#F28C28\", \"Endothelial\" = \"#F28C28\", \"Muscle\" = \"#DD3F4E\", \"Mesenchymal\" = \"#FFB6C1\", \"Mesenchymal\" = \"#FFB6C1\", \"Mesenchymal\" = \"#FFB6C1\", \"Mesenchymal\" = \"#FFB6C1\", \"Mesenchymal\" = \"#FFB6C1\", \"Mesenchymal\" = \"#FFB6C1\", \"Mesenchymal\" =\"#FFB6C1\")\n",
    "cell_counts <- as.numeric(table(factor(combined$cluster_anno_l2, levels = cell_order)))\n",
    "names(cell_order) <- cal2_cols\n",
    "anno_df <- data.frame(cell_order, cell_lineages)\n",
    "names(cal2_cols) <- levels(droplevels(as.factor(combined$cluster_anno_l2)))\n",
    "ha = HeatmapAnnotation(annotation_label = c(\"Cell Type\", \"Cell Lineage\"), df = anno_df,\n",
    "                       border = TRUE, which = 'col', col = list(cell_order = cal2_cols,cell_lineages = c(cell_lineage_colors)))\n",
    "avg <- as.data.frame(avg$RNA)\n",
    "genes_show = c(\"AVP\",\"SPINK2\",\"IL7R\", \"GZMA\",\"DNTT\",\"CD79B\",\"PAX5\",\"MZB1\", \"MPO\", \"SRGN\",\"S100A9\" ,\"LYZ\", \"CD14\", \"C1Q\",\"GATA1\", \"ITGA2B\", \"HBA1\",\"HBA2\", \"CXCL12\",\"LEPR\", \"PDGFRA\", \"NCAM1\",\"FLT1\",\"EMCN\", \"TAGLN\", \"ACTA2\")\n",
    "Heatmap(t(scale(t(avg[,cell_order]))), name = \"mat\", rect_gp = gpar(col = \"black\", lwd = 0), col = col_fun, column_order = cell_order,\n",
    "        column_title = \"scRNA Average Scaled Cell Type Marker Expression\", clustering_method_rows = \"single\", show_row_names = FALSE, row_names_gp = grid::gpar(fontsize = 8), top_annotation =  ha) +\n",
    "  rowAnnotation(link=anno_mark(at=which(rownames(avg) %in% genes_show), labels = rownames(avg)[rownames(avg) %in% genes_show], which = \"rows\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ab59475-bc64-4eb3-9a2b-4a21b49ad1a8",
   "metadata": {},
   "source": [
    "# Perform Analysis of MSC and Endothelial Cell Subsets -----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eadd6422-9c08-4352-9c00-663118f4d1e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subset MSCs, Re-Normalize/Scale, Compute PCA and UMAP, leaving out the RNAlo MSCs ----\n",
    "MSCs <- subset(combined, cluster_anno_l2 == \"Fibro-MSC\" | cluster_anno_l2 == \"THY1+ MSC\" | cluster_anno_l2 == \"Adipo-MSC\" | cluster_anno_l2 == \"OsteoFibro-MSC\" | cluster_anno_l2 == \"Osteo-MSC\" | cluster_anno_l2 == \"Osteoblast\")\n",
    "MSCs <- NormalizeData(MSCs)\n",
    "MSCs <- FindVariableFeatures(MSCs)\n",
    "MSCs <- ScaleData(MSCs)\n",
    "MSCs <- RunPCA(MSCs, reduction.name = \"MSC_pca\", features = VariableFeatures(MSCs))\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:10, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim10\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:20, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim20\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:30, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim30\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:40, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim40\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:50, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim50\")\n",
    "\n",
    "DimPlot(MSCs, label = TRUE, reduction = \"MSC_UMAP_dim50\") + coord_fixed() + NoAxes() + NoLegend()\n",
    "VlnPlot(subset(combined, cluster_anno_l2 == \"Fibro-MSC\" | cluster_anno_l2 == \"THY1+ MSC\" | cluster_anno_l2 == \"Adipo-MSC\" | cluster_anno_l2 == \"OsteoFibro-MSC\" | cluster_anno_l2 == \"Osteo-MSC\" | cluster_anno_l2 == \"Osteoblast\"\n",
    "               | cluster_anno_l2 == \"RNAlo MSC\"), features = c(\"nCount_RNA\", 'nFeature_RNA', 'percent.mt'), pt.size = 0, group.by = 'cluster_anno_l2')\n",
    "\n",
    "MSCs$original_seurat_clusters <- MSCs@active.ident\n",
    "MSCs <- FindNeighbors(MSCs, dims = 1:20)\n",
    "MSCs <- FindClusters(MSCs, algorithm = 2, resolution = 0.25, group.singletons = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afcda1f7-89d5-4519-9ba3-415b1f74628b",
   "metadata": {},
   "source": [
    "# update viscello with annotated combined obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77b33ec8-1b7b-47a7-9af6-774c55be60c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "DefaultAssay(combined) <- \"RNA\"\n",
    "fmeta <- data.frame(symbol = rownames(combined)) \n",
    "rownames(fmeta) <- fmeta$symbol\n",
    "eset <- new(\"ExpressionSet\",\n",
    "            assayData = assayDataNew(\"environment\", exprs=combined@assays$RNA@counts, \n",
    "                                     norm_exprs = combined@assays$RNA@data),\n",
    "            phenoData =  new(\"AnnotatedDataFrame\", data = combined@meta.data),\n",
    "            featureData = new(\"AnnotatedDataFrame\", data = fmeta))\n",
    "saveRDS(eset, \"VisCello/eset.rds\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04bc0fd6-7f5e-4489-9d07-d41959bc9527",
   "metadata": {},
   "source": [
    "# Creating a cello for all the cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef1cfab0-90cb-4aef-b6da-3d95cc5d522b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cello <- new(\"Cello\", name = \"Combined All Cells SB66\", idx = 1:ncol(eset)) # Index is basically the column index, here all cells are included \n",
    "# Code for computing dimension reduction, not all of them is necessary, and you can input your own dimension reduction result into the cello@proj list.\n",
    "# It is also recommended that you first filter your matrix to remove low expression genes and cells, and input a matrix with variably expressed genes\n",
    "cello@proj <- list('PCA' = combined@reductions$pca@cell.embeddings,\n",
    "                   'UMAP_dim30' =combined@reductions$UMAP_dim30@cell.embeddings,\n",
    "                   'UMAP_dim50' =combined@reductions$UMAP_dim50@cell.embeddings)\n",
    "clist <- list()\n",
    "clist[[\"SB66 Global dataset\"]] <- cello\n",
    "saveRDS(clist, \"VisCello/clist.rds\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "876f9364-8340-4c9d-ad6f-321cb3f37255",
   "metadata": {},
   "source": [
    "# Add zoom of MSC to Viscello "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f628678e-6e53-4822-9a62-87704d9278fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "zoom_type <- \"MSC_Trajectory\"\n",
    "cur_idx <-  which(colnames(eset) %in% colnames(MSCs))\n",
    "cello <- new(\"Cello\", name = zoom_type, idx = cur_idx) \n",
    "cello@proj <- list('MSC_PCA' = MSCs@reductions$MSC_pca@cell.embeddings,\n",
    "                   'MSC_UMAP_dim10'=MSCs@reductions$MSC_UMAP_dim10@cell.embeddings, \n",
    "                   'MSC_UMAP_dim20'=MSCs@reductions$MSC_UMAP_dim20@cell.embeddings, \n",
    "                   'MSC_UMAP_dim30'=MSCs@reductions$MSC_UMAP_dim30@cell.embeddings, \n",
    "                   'MSC_UMAP_dim40'=MSCs@reductions$MSC_UMAP_dim40@cell.embeddings,\n",
    "                   'MSC_UMAP_dim50'=MSCs@reductions$MSC_UMAP_dim50@cell.embeddings)\n",
    "clist[[zoom_type]] <- cello\n",
    "saveRDS(clist, \"VisCello/clist.rds\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aa45e55-07d1-4974-bd23-da221d235ee3",
   "metadata": {},
   "source": [
    "# remove doublets (PTPRC+ cluster) identified in viscello and recompute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bb19e46-6e41-4f5e-830e-2442387fd78e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "doublet_ids <- read_csv(\"MSC_Doublets_Data_VisCello.csv\")\n",
    "MSCs <- subset(MSCs, invert = TRUE, cells = doublet_ids$...1)\n",
    "# recompute umap\n",
    "MSCs <- NormalizeData(MSCs)\n",
    "MSCs <- FindVariableFeatures(MSCs)\n",
    "MSCs <- ScaleData(MSCs)\n",
    "MSCs <- RunPCA(MSCs, reduction.name = \"MSC_pca\", features = VariableFeatures(MSCs))\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:10, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim10\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:20, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim20\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:30, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim30\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:40, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim40\")\n",
    "MSCs <- RunUMAP(MSCs, dims = 1:50, reduction = \"MSC_pca\", reduction.name = \"MSC_UMAP_dim50\")\n",
    "MSCs <- FindNeighbors(MSCs, dims = 1:20)\n",
    "MSCs <- FindClusters(MSCs, algorithm = 2, resolution = 0.25, group.singletons = FALSE)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4aaba99-77ac-42e3-abde-2b00d18e0a1d",
   "metadata": {},
   "source": [
    "# Repeat for Endothelial Cells ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c6ba6dc-ad01-4519-b3b3-51979210b1d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subset Endothelial Cells, Re-Normalize/Scale, Compute PCA and UMAP ----\n",
    "Endo <- subset(combined, cluster_anno_l2 == \"AEC\" | cluster_anno_l2 == \"SEC\")\n",
    "Endo <- NormalizeData(Endo)\n",
    "Endo <- FindVariableFeatures(Endo)\n",
    "Endo <- ScaleData(Endo)\n",
    "Endo <- RunPCA(Endo, reduction.name = \"Endo_pca\", features = VariableFeatures(Endo))\n",
    "Endo <- RunUMAP(Endo, dims = 1:10, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim10\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:20, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim20\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:30, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim30\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:40, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim40\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:50, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim50\")\n",
    "Endo <- FindNeighbors(Endo, dims = 1:20)\n",
    "Endo <- FindClusters(Endo, algorithm = 2, resolution = 0.25, group.singletons = FALSE)\n",
    "Endo_Markers <- FindAllMarkers(Endo)\n",
    "write.csv(Endo_Markers, \"Endo_Markers.csv\")\n",
    "\n",
    "DimPlot(Endo, label = TRUE, reduction = \"Endo_UMAP_dim50\") + coord_fixed() + NoAxes() + NoLegend()\n",
    "VlnPlot(Endo, features = c(\"nCount_RNA\", 'nFeature_RNA', 'percent.mt'), pt.size = 0, group.by = 'seurat_clusters') # note high nCount and nFeature for clusters 4 and 5, extremely low for Cluster 3\n",
    "FeaturePlot(Endo, reduction = \"Endo_UMAP_dim50\", features= c(\"CDH5\", \"PTPRC\", \"CD34\", \"DCN\", \"CXCL12\", \"KITLG\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83dc763c-bd31-488c-855f-77f90ff3d982",
   "metadata": {},
   "source": [
    "# Add the new zoom of Endo with putative doublets removed to Viscello"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a30e7b42-81ad-46ec-af7b-e5d3be0d5215",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "zoom_type <- \"Endo_Trajectory\"\n",
    "cur_idx <-  which(colnames(eset) %in% colnames(Endo))\n",
    "cello <- new(\"Cello\", name = zoom_type, idx = cur_idx) \n",
    "cello@proj <- list('Endo_PCA' = Endo@reductions$Endo_pca@cell.embeddings,\n",
    "                   'Endo_UMAP_dim10'=Endo@reductions$Endo_UMAP_dim10@cell.embeddings, \n",
    "                   'Endo_UMAP_dim20'=Endo@reductions$Endo_UMAP_dim20@cell.embeddings, \n",
    "                   'Endo_UMAP_dim30'=Endo@reductions$Endo_UMAP_dim30@cell.embeddings, \n",
    "                   'Endo_UMAP_dim40'=Endo@reductions$Endo_UMAP_dim40@cell.embeddings,\n",
    "                   'Endo_UMAP_dim50'=Endo@reductions$Endo_UMAP_dim50@cell.embeddings)\n",
    "clist[[zoom_type]] <- cello\n",
    "saveRDS(clist, \"VisCello/clist.rds\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d4fc44-7b4a-4094-98ec-d880ce8b7e57",
   "metadata": {},
   "source": [
    "# remove doublets (PTPRC+ and DCN+ cluster) identified in viscello and recompute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5120e64f-d3f1-46fe-a1d5-f5041b9cbc99",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Endo_doublet_ids <- read_csv(\"~/Documents/NBM_Microenvironment/NBM_Atlas_scRNA/Final_scRNA_Analysis/Endo_Doublets_LQcells_Removed_Data_VisCello.csv\")\n",
    "Endo <- subset(Endo,  cells = Endo_doublet_ids$...1) # remove doublets\n",
    "Endo <- subset(Endo, subset = seurat_clusters != 3) # remove likely low quality cells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fefe61d6-e80d-4473-af00-cd0587bc54e1",
   "metadata": {},
   "source": [
    "# recompute umap and clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d40fa3cd-58b6-415a-afc5-fcd3ada6d404",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Endo <- NormalizeData(Endo)\n",
    "Endo <- FindVariableFeatures(Endo)\n",
    "Endo <- ScaleData(Endo)\n",
    "Endo <- RunPCA(Endo, reduction.name = \"Endo_pca\", features = VariableFeatures(Endo))\n",
    "Endo <- RunUMAP(Endo, dims = 1:10, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim10\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:20, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim20\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:30, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim30\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:40, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim40\")\n",
    "Endo <- RunUMAP(Endo, dims = 1:50, reduction = \"Endo_pca\", reduction.name = \"Endo_UMAP_dim50\")\n",
    "Endo <- FindNeighbors(Endo, dims = 1:20)\n",
    "Endo <- FindClusters(Endo, algorithm = 2, resolution = 1, group.singletons = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90dfd3aa-ebc5-44c7-bb9d-c64a7f2a91e7",
   "metadata": {},
   "source": [
    "# Figure 2A Make UMAP for MSC Subsets ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1349169a-518d-4b1b-8556-41395abcfe67",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p1 <- DimPlot(MSCs, label = TRUE, group.by = 'cluster_anno_l2',reduction = \"MSC_UMAP_dim50\", raster = TRUE, raster.dpi = c(4112,4112), pt.size = 8, cols = cal2_cols) + coord_fixed() + NoAxes() + NoLegend()\n",
    "p1_raster <- rasterize(p1)\n",
    "ggsave(p1_raster,filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Figure2/panel/PanelA_scRNA_MSC_UMAP.eps\", device = 'eps', width = 5,height = 5)\n",
    "ggsave(p1_raster,filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Figure2/panel/PanelA_scRNA_MSC_UMAP.pdf\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "972558dd-e11c-4396-ba2c-e7cc39a07400",
   "metadata": {},
   "source": [
    "## Supplemental Figure S3A RNAlo MSC QC ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0418ec2-6b5b-4e8a-9571-3f162206fb13",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "VlnPlot(subset(combined, cluster_anno_l2 %in% c(\"Adipo-MSC\", \"Osteo-MSC\", \"THY1+ MSC\", \"OsteoFibro-MSC\", \"Osteoblast\", \"Fibro-MSC\", \"RNAlo MSC\")), features = c(\"nFeature_RNA\", \"nCount_RNA\", \"percent.mt\"), pt.size = 0, cols = cal2_cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b947533-fc58-4b4c-be81-ef40eb6ec49e",
   "metadata": {},
   "source": [
    "## Supplemental Figure S3B - Various Canonical MSC genes ---- "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62cec4ef-6b4d-46e0-b55f-9ecd2f1d5c20",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "mesenchymal_cell_types <-  c(\"Adipo-MSC\",\"THY1+ MSC\",\"Fibro-MSC\",\"OsteoFibro-MSC\",\"Osteo-MSC\",\"Osteoblast\", \"SEC\", \"AEC\", \"VSMC\") # need to reload the original untransformed combined object\n",
    "mes <- subset(combined, cluster_anno_l2 %in% mesenchymal_cell_types)\n",
    "mes$cluster_anno_l2 <- factor(mes$cluster_anno_l2, levels =c(\"Fibro-MSC\",\"OsteoFibro-MSC\",\"Osteo-MSC\",\"Osteoblast\",\"Adipo-MSC\",\"THY1+ MSC\",\"AEC\", \"SEC\", \"VSMC\"))\n",
    "p1 <- DotPlot(mes,group.by = \"cluster_anno_l2\", scale = FALSE, features = c(\"NES\",\"LBP\",\"NCAM1\", \"SPP1\", \"DCN\", \"CSPG4\"), cluster.idents = FALSE) & RotatedAxis() & scale_color_viridis(option = \"plasma\") & coord_flip()\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelS2B_DotPlot_ExtraMarkerGenes_v1.pdf', width = 8, height = 3.5)\n",
    "\n",
    "\n",
    "\n",
    "DotPlot(subset(combined, cluster_anno_l2 %in% mesenchymal_cell_types),scale = FALSE, cols = \"Inferno\", features = c(\"NES\",\"LBP\",\"NCAM1\", \"SPP1\", \"DCN\", \"CSPG4\")) & RotatedAxis() \n",
    "p1 <- VlnPlot(subset(combined, cluster_anno_l2 %in% mesenchymal_cell_types),slot = 'data' , pt.size = 0,  features = \"NES\", cols = cal2_cols) & RotatedAxis() & NoLegend()\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelB_NES_VlnPlots.pdf',\n",
    "       device ="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "226b7ce5-2f2b-4211-a68e-38b2b2f64cd3",
   "metadata": {},
   "source": [
    "## Supplemental Figure S3C MSC and Endo frequency per sample ----#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c313df8-2485-4673-8a2e-364e125b8bcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "MSCs$cluster_anno_l2 <- droplevels(MSCs$cluster_anno_l2)\n",
    "cell_counts_MSC <- as.data.frame(table(MSCs$cluster_anno_l2,MSCs$orig.ident))\n",
    "\n",
    "cell_counts_MSC$Var2 <- droplevels(cell_counts_MSC$Var1)\n",
    "MSC_cols <- MetBrewer::met.brewer(\"Austria\", n = 6)\n",
    "MSC_cols <- RColorBrewer::brewer.pal(n = 6, name = \"Set2\")\n",
    "\n",
    "\n",
    "p3 <- ggplot(data=cell_counts_MSC, aes(x=forcats::fct_rev(Var2),y=Freq, fill=Var1)) +\n",
    "  geom_bar(position=\"fill\",stat=\"identity\") + \n",
    "  coord_flip() +\n",
    "  theme_bw() + \n",
    "  labs(y='Cell Type Frequency') + \n",
    "  labs(x='') + scale_fill_manual(values=MSC_cols) + theme(axis.text = element_text(size=16)) + ggtitle(\"CODEX Cell Type Frequencies Per Sample\") + scale_x_discrete(labels = c(\"H41\", \"H39\", \"H38\", \"H36\", \"H35\", \"H34\", \"H33\", \"H32\", \"H24\", \"H23\", \"H21\", \"H14\") )\n",
    "\n",
    "Endo$cluster_anno_l2 <- droplevels(Endo$cluster_anno_l2)\n",
    "cell_counts_Endo <- as.data.frame(table(Endo$cluster_anno_l2,Endo$orig.ident))\n",
    "\n",
    "\n",
    "\n",
    "p4 <- ggplot(data=cell_counts_Endo, aes(x=forcats::fct_rev(Var2),y=Freq, fill=Var1)) +\n",
    "  geom_bar(position=\"fill\",stat=\"identity\") + \n",
    "  coord_flip() +\n",
    "  theme_bw() + \n",
    "  labs(y='Cell Type Frequency') + \n",
    "  labs(x='') + scale_fill_manual(values= c(\"#AA336A\", \"#93E9BE\")) + theme(axis.text = element_text(size=16)) + ggtitle(\"CODEX Cell Type Frequencies Per Sample\") + scale_x_discrete(labels = c(\"H41\", \"H39\", \"H38\", \"H36\", \"H35\", \"H34\", \"H33\", \"H32\", \"H24\", \"H23\", \"H21\", \"H14\") + RotatedAxis())\n",
    "\n",
    "ggsave(p3, device = \"pdf\", height = 3, width = 4, units = \"in\", filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/MSC_Frequency.pdf\")\n",
    "ggsave(p4, device = \"pdf\", height = 3, width = 4, units = \"in\", filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/Endo_Frequency.pdf\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9fa7944-1f95-4c61-92a5-fde8555a20fe",
   "metadata": {},
   "source": [
    "## calculate relative frequency of each cell type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "feb8bc33-436e-477c-9dcd-2114d482892a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cell_counts_MSC %>% dplyr::group_by(Var1) %>% summarise(rel_freq = Freq/n())\n",
    "\n",
    "cell_counts_MSC <- cell_counts_MSC %>%\n",
    "  group_by(Var2) %>%\n",
    "  mutate(total = sum(Freq))\n",
    "\n",
    "cell_counts_MSC$rel_freq <-  cell_counts_MSC$Freq / cell_counts_MSC$total\n",
    "cell_counts_MSC_summarystats <-  cell_counts_MSC %>% group_by(Var1) %>% summarise(median=median(rel_freq), std.dev = sd(rel_freq), min = min(rel_freq), max = max(rel_freq))\n",
    "cell_counts_MSC %>% group_by(Var2) %>% summarise(sum(rel_freq))\n",
    "\n",
    "cell_counts_MSC$is_thy1oradipo <- cell_counts_MSC$Var1 == \"THY1+ MSC\" | cell_counts_MSC$Var1 == \"Adipo-MSC\"\n",
    "cell_counts_MSC_summarystats <-  cell_counts_MSC %>% group_by(is_thy1oradipo) %>% summarise(median=median(rel_freq), std.dev = sd(rel_freq), min = min(rel_freq), max = max(rel_freq))\n",
    "\n",
    "cell_counts_Endo <- cell_counts_Endo %>%\n",
    "  group_by(Var2) %>%\n",
    "  mutate(total = sum(Freq))\n",
    "\n",
    "cell_counts_Endo$rel_freq <-  cell_counts_Endo$Freq / cell_counts_Endo$total\n",
    "cell_counts_Endo_summarystats <-  cell_counts_Endo %>% group_by(Var1) %>% summarise(median=median(rel_freq), std.dev = sd(rel_freq), min = min(rel_freq), max = max(rel_freq))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ddbaeb-670e-4b0d-9b31-e4ebf5fe0f2d",
   "metadata": {},
   "source": [
    "# test that SEC is more than AEC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae787775-1279-4316-83b4-1ea312585b74",
   "metadata": {},
   "outputs": [],
   "source": [
    "prop.test(7,12, p=0.5, correct = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36977164-8ce7-4db7-9896-ed84b2b5abcb",
   "metadata": {},
   "source": [
    "## Supplemental Figure S4A - CSF3 -----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "699fbd54-478d-4a1b-be9a-46fec4f0a85e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p1 <- VlnPlot(combined, features = c(\"CSF3\"), pt.size = 1,  sort = TRUE) + NoLegend()\n",
    "ggsave(p1, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelD_GCSF_markers.pdf\", device = \"pdf\", height = 5, width = 12)\n",
    "\n",
    "p1 <- FeaturePlot(combined, max.cutoff = 'q99', features = c(\"CSF3\"), coord.fixed = TRUE, order = TRUE, cols=brewer.pal(n = 100, name = \"Reds\")) & NoAxes() # note use of order=TRUE to highlight the positive cells\n",
    "p1_raster <- rasterize(p1, dpi = 300)\n",
    "ggsave(p1_raster, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelD_CSF3_UMAP.pdf\", device = \"pdf\", height = 5, width = 7)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77df599a-19ff-40c8-9c41-588e61dc5764",
   "metadata": {},
   "source": [
    "# Not a figure - but sorting strategy justification, related to Fig S3E -----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3653e98f-e546-4445-8544-e35f4d58ba3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "DotPlot(MSCs, group.by = 'cluster_anno_l2', features = c(\"PTPRC\",\"GYPA\", \"CD38\",\"CDH5\", \"PDPN\", \"NCAM1\", \"LEPR\", \"THY1\"), scale = FALSE, cluster.idents = TRUE) + scale_color_viridis(option = 'plasma') + RotatedAxis()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b13d963-fc48-451c-b880-5c04f72f485c",
   "metadata": {},
   "source": [
    "# Supplemental Table 2 - calculate DEGs between MSC and Endo subsets ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ba3696a-cc0a-42b8-99e5-8fa7d87fccc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "MSCs<- SetIdent(MSCs, value = \"cluster_anno_l2\")\n",
    "MSC_markers <- FindAllMarkers(MSCs)\n",
    "Endo_markers <- FindAllMarkers(Endo)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21ebff5d-90a5-4111-bc1c-022368d7b08e",
   "metadata": {},
   "source": [
    "# Figure 2B - Dot Plot of canonical mesenchymal genes ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "996f5d4a-4049-4c7c-bbc9-7b9d1a6ad331",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 <- DotPlot(MSCs, group.by = 'cluster_anno_l2', features = c(\"PDGFRA\", \"VIM\",\"CXCL12\", \"LEPR\", \"LPL\", \"APOE\",\"THY1\",\"RUNX2\", \"SP7\", \"IBSP\", \"BGLAP\", \"COL1A1\",\"GSN\", \"APOD\", \"PDPN\", \"HAS1\",\"DPT\", \"CD164\", \"NT5E\"), scale = FALSE, cluster.idents = FALSE) + scale_color_viridis(option = 'plasma') + RotatedAxis()\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Figure2/panel/PanelB_DotPlot_MCS_DEGs.eps',\n",
    "       device = 'pdf', width = 10, height = 3.5)\n",
    "# Figure 2C Dot Plot of ISCT genes and commonly used MSC markers ----\n",
    "p1 <- DotPlot(MSCs, group.by = 'cluster_anno_l2', features = c(\"NGFR\", \"MCAM\", \"NT5E\", \"THY1\", \"ENG\"), scale = FALSE, cluster.idents = FALSE) + scale_color_viridis(option = 'plasma') + RotatedAxis()\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Figure2/panel/PanelB_DotPlot_MCS_DEGs.eps',\n",
    "       device = 'pdf', width = 10, height = 3.5)\n",
    "\n",
    "# Figure 2H Make dot plot of DEGs between SEC and AEC -----\n",
    "Endo_markers <- FindMarkers(Endo, group.by = 'cluster_anno_l2', ident.1 = \"AEC\", ident.2 = \"SEC\")\n",
    "Endo_markers$gene_name <- rownames(Endo_markers)\n",
    "Endo_markers%>% \n",
    "  top_n(n = 10, wt = avg_log2FC) -> top10_AEC # Get 10 most significant genes for each cell type\n",
    "Endo_markers  %>% \n",
    "  top_n(n = 10, wt = -avg_log2FC) -> top10_SEC # Get 10 most significant genes for each cell type\n",
    "# note chose both some canonical markers and top DEGs\n",
    "p1 <- DotPlot(Endo, group.by = 'cluster_anno_l2', features = c(\"CDH5\", \"CD34\", \"KDR\",\"PODXL\", \"ICAM2\", rownames(top10_AEC), \"NR2F2\", \"EPHB4\", rownames(top10_SEC)), scale = FALSE,  cluster.idents = TRUE) + scale_color_viridis(option = 'plasma') + RotatedAxis()\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Figure2/panel/PanelD_DotPlot_BMEC_DEGs.pdf', width = 10, height = 3.5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91f11946-cba3-4136-8bab-2383d678516d",
   "metadata": {},
   "source": [
    "# Figure 3A - Plot Notable Putative Hematopoietic Supportive Factors ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f265a1-12ad-43a8-866b-5cbf2807884e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cytokines <- c(\"TGFB1\",\"TNFRSF11B\",\"CDH2\",\"SPP1\", \"DLL1\", \"DLL4\", \"JAG1\", \"JAG2\",\"SELE\", \"SELP\",\"SELPLG\",\"FLT3LG\", \"CSF1\", \"CXCL12\", \"KITLG\",\"VCAM1\", \"TNFSF11\", \"IL6\",\"IL7\",\"PTN\", \"FGF1\", \"FGF2\",\"NGF\", \"ANGPT1\", \"ANGPT2\", \"IGF1\")\n",
    "mesenchymal_cell_types <-  c(\"Adipo-MSC\",\"THY1+ MSC\",\"Fibro-MSC\",\"OsteoFibro-MSC\",\"Osteo-MSC\",\"Osteoblast\", \"SEC\", \"AEC\", \"VSMC\") # need to reload the original untransformed combined object\n",
    "mes <- subset(combined, cluster_anno_l2 %in% mesenchymal_cell_types)\n",
    "mes$cluster_anno_l2 <- factor(mes$cluster_anno_l2, levels =c(\"Fibro-MSC\",\"OsteoFibro-MSC\",\"Osteo-MSC\",\"Osteoblast\",\"Adipo-MSC\",\"THY1+ MSC\",\"AEC\", \"SEC\", \"VSMC\"))\n",
    "p1 <- DotPlot(mes,group.by = \"cluster_anno_l2\", scale = FALSE,  features = cytokines, cluster.idents = FALSE) & RotatedAxis() & scale_color_viridis(option = \"plasma\")\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Figure3/panels/PanelA_DotPlot_SupportiveFactors_v3.pdf', width = 10, height = 3.5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b6e3e9c-cb8d-4f29-8bbb-70b7db70f192",
   "metadata": {},
   "source": [
    "# Figure S4B Fetal MSC Reference Mapping ------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff478af4-56b3-47ee-98ff-32b1a09e68bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# With fetal bone marrow\n",
    "Convert(\"fig1b_fbm_scaled_gex_updated_dr_20210104.h5ad\", dest = \"h5seurat\", overwrite = TRUE) # Convert h5ad obj from Jardine paper to seurat\n",
    "fbm <- LoadH5Seurat(\"~/Documents/NBM_Microenvironment/NBM_Atlas_scRNA/Final_scRNA_Analysis/Objects_To_RefMap/fig1b_fbm_scaled_gex_updated_dr_20210104.h5seurat\")\n",
    "fbm_stroma <- subset(fbm, broad_fig1_cell.labels == \"stroma\")\n",
    "fbm_stroma_MSC <- subset(fbm_stroma, cell.labels == \"adipo-CAR\" | cell.labels == \"osteoblast\" | cell.labels == \"osteoblast precursor\" | cell.labels == \"osteochondral precursor\" | cell.labels ==\"endosteal fibroblast\" | cell.labels ==\"arteriolar fibroblast\") # chondrocytes and muscle/muscle stem cells not included because they aren't represented in our dataset (not in adult BM)\n",
    "\n",
    "remove(fbm)\n",
    "anchors <- FindTransferAnchors(\n",
    "  reference = MSCs,\n",
    "  query = fbm_stroma_MSC,\n",
    "  normalization.method = \"LogNormalize\",\n",
    "  reference.reduction = \"MSC_pca\",\n",
    "  dims = 1:50\n",
    ")\n",
    "fbm_stroma_MSC <- MapQuery(\n",
    "  anchorset = anchors,\n",
    "  query = fbm_stroma_MSC,\n",
    "  reference = MSCs,\n",
    "  refdata = list(\n",
    "    MSC_refmap = \"cluster_anno_l2\"\n",
    "  ),\n",
    "  reference.reduction = \"MSC_pca\", \n",
    "  reduction.model = \"MSC_UMAP_dim50\"\n",
    ")\n",
    "\n",
    "DimPlot(fbm_stroma_MSC, reduction = \"umap\", group.by = \"predicted.MSC_refmap\", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()\n",
    "p2 <- ggplot(fbm_stroma_MSC@meta.data, aes(fill=predicted.MSC_refmap, x=orig.ident)) + \n",
    "  geom_bar(position=\"fill\", stat=\"count\") + theme_minimal() + RotatedAxis() +\n",
    "  ylab(\"Percentage of Total MSCs\") + xlab(\"Sample\") + scale_fill_manual(values=cal2_cols, limits=force)\n",
    "\n",
    "plot_grid(p1+p2)\n",
    "ggsave(p1, filename = '~/Documents/Manuscripts/NBM_Atlas/Figures/Figure2/panel/PanelG_StackedBar_Fetal_MSCs.pdf',width = 5, height = 5)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "442d6e77-287b-44ba-a591-3a3799adf821",
   "metadata": {},
   "source": [
    "## Supplemental Figure S4E - Lymphatic Endothelial cell markers ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b919c853-ceaf-4842-b384-15f32835af9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p1 <- FeaturePlot(Endo, reduction = \"Endo_UMAP_dim30\", features = c(\"LYVE1\"), pt.size = 0.1, max.cutoff = 'q99', cols=brewer.pal(n = 100, name = \"Reds\"), coord.fixed = TRUE, ncol = 3) & NoAxes()\n",
    "p1_raster <- rasterize(p1, dpi = 300)\n",
    "ggsave(p1_raster, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelC_LEC_markers_LYVE1.pdf\", device = \"pdf\", height = 5, width = 7)\n",
    "p2 <- FeaturePlot(Endo, reduction = \"Endo_UMAP_dim30\", features = c(\"PROX1\"), pt.size = 0.1, max.cutoff = 'q99', cols=brewer.pal(n = 100, name = \"Reds\"), coord.fixed = TRUE, ncol = 3) & NoAxes()\n",
    "p2_raster <- rasterize(p2, dpi = 300)\n",
    "ggsave(p2_raster, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelC_LEC_markers_PROX1.pdf\", device = \"pdf\", height = 5, width = 7)\n",
    "p3 <- FeaturePlot(Endo, reduction = \"Endo_UMAP_dim30\", features = c(\"PDPN\"), pt.size = 0.1, max.cutoff = 'q99', cols=brewer.pal(n = 100, name = \"Reds\"), coord.fixed = TRUE, ncol = 3) & NoAxes()\n",
    "p3_raster <- rasterize(p3, dpi = 300)\n",
    "ggsave(p3_raster, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure2/PanelC_LEC_markers_PDPN.pdf\", device = \"pdf\", height = 5, width = 7)\n",
    "\n",
    " \n",
    "FeaturePlot(MSCs, reduction = \"MSC_UMAP_dim30\", features = c(\"LYVE1\", \"PROX1\", \"PDPN\"), pt.size = 0.1, max.cutoff = 'q99', coord.fixed = TRUE, cols=brewer.pal(n = 100, name = \"Reds\")) & NoAxes()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d944f82-ae4a-4dbb-b86e-35c848f520e1",
   "metadata": {},
   "source": [
    "# Figure 5E - AUCell Hypoxia Score Analysis ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a8af947-7f8f-4eb6-8d21-e71a5753cc28",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "myeloid_development <- c(\"HSC\", \"MPP\",\"Cycling HSPC\", \"GMP\", \"Early Myeloid Progenitor\", \"Late Myeloid\", \"Neutrophil\")\n",
    "FeaturePlot(combined_AUC, reduction = \"UMAP_dim30\", features = \"HALLMARK_HYPOXIA\", coord.fixed = TRUE, max.cutoff = 'q99') + \n",
    "  scale_color_gradientn(colours = rev(brewer.pal(n = 11, name = \"RdBu\"))) + NoAxes()\n",
    "\n",
    "VlnPlot(subset(combined_AUC, subset = cluster_anno_l2 %in% myeloid_development), group.by = 'cluster_anno_l2', features = \"HALLMARK_HYPOXIA\", pt.size = 0, sort = TRUE, cols = cal2_cols) + NoLegend()\n",
    "\n",
    "combined_AUC_myeloid <- combined_AUC %>% subset(cluster_anno_l2 %in% myeloid_development)\n",
    "\n",
    "# Perform one vs. rest T.test with multiple hypothesis correction for each population\n",
    "test <- c()\n",
    "c <- c()\n",
    "d <- c()\n",
    "for (i in levels(droplevels(combined_AUC_myeloid$cluster_anno_l2))) {\n",
    "  a <- subset(combined_AUC_myeloid@meta.data, cluster_anno_l2 == {i})$HALLMARK_HYPOXIA\n",
    "  d <- append(x = d, values = i)\n",
    "  b <- subset(combined_AUC_myeloid@meta.data, cluster_anno_l2 != {i})$HALLMARK_HYPOXIA\n",
    "  test <- t.test(a,b, alternative = \"less\")\n",
    "  c <- append(x = c, values = test$p.value)\n",
    "}\n",
    "#correct for multiple hypothesis \n",
    "c <- p.adjust(c,method = \"BH\")\n",
    "names(c) <- d\n",
    "c\n",
    "\n",
    "subset(combined_AUC_myeloid@meta.data, cluster_anno_l2 == \"Neutrophil\")$HALLMARK_HYPOXIA\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16ad6363-9a1e-4e43-8858-6021f0491b58",
   "metadata": {},
   "source": [
    "# Supplemental Figure S7B Distance to Bone ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a28b24e-2e79-476a-9946-a75356b0a4f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "cn_ranks <- read_csv(\"~/Documents/NBM_Microenvironment/NBM_Atlas_CODEX_Samples/Non-Cell Microenvironment Analysis/combined_neighbor.csv\")\n",
    "\n",
    "nbs_names <- c(\"HSC / Mature Myeloid\", \"Erythroid/Myeloid\", \"PC/Arteriolar\", \"Erythroid\", \"Arteriolar\", \"Erythroid\", \"Lymphoid\", \"Erythroid/Myeloid/Lymphoid\", \"Early Myeloid / Endosteal\", \"Myeloid/Lymphoid\", \"HSPC/Intermediate Myeloid\", \"Erythroid/Myeloid/Lymphoid\", \"Erythroid/Myeloid\", \"Early Myeloid / Arteriolar\", \"Peri-Arterolar Lymphoid\")\n",
    "nbs_cols <- c(\"#A8D37F\", \"#51C8EB\", \"#FF00FF\", \"#FBED24\", \"#FF0000\", \"#FBED24\", \"#89919C\", \"#BD7CB5\", \"#4B409A\", \"#FFA500\", \"#A15A26\", \"#BD7CB5\", \"#51C8EB\", \"#4B409A\", \"#FF007F\")\n",
    "names(nbs_cols) <- nbs_names\n",
    "\n",
    "p2 <- cn_ranks %>% dplyr::filter(cn_ranks$structure == \"bone\") %>% ggplot(aes(x = reorder(neighborhoods, -normalized_rank), y = normalized_rank, fill = neighborhoods)) + geom_boxplot() + theme_minimal() + NoLegend()+ RotatedAxis() + scale_fill_manual(values = nbs_cols)\n",
    "ggsave(p2, filename = \"~/Documents/Manuscripts/NBM_Atlas/Figures/Supplemental_Figures/Related_to_Figure5/neighborhoods_bonedist_normrank.pdf\", device = \"pdf\", units = \"in\", width = 8, height = 4)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R(conda R_sing11)",
   "language": "R",
   "name": "r_sing1"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
